import json

import laspy
from rasterio.transform import from_origin
import open3d as o3d
from tools.LoggerTools import *
import cv2
from PIL import ImageDraw, ImageFont

import rasterio
from rasterio.plot import reshape_as_image
from rasterio.transform import xy
from PIL import Image, ImageDraw
import numpy as np
import matplotlib.pyplot as plt


setup_global_logging()


# 调整参数
denose_nb_neighbors = 80  # 噪声的相邻要素数
denose_std_ratio = 0.5     # 噪声去除的标准偏差比率
resolution = 0.5           # DEM分辨率，值越小越精细但慢（例如 0.5=每个像素显示0.5平方米的高程值）
threshold = 0.2            # 高程变化阈值（以米为单位）


def get_file_paths(base_path):
    if not os.path.exists(base_path):
        raise FileNotFoundError(f"路径不存在: {base_path}")
        # 确保是目录
    if not os.path.isdir(base_path):
        raise NotADirectoryError(f"路径不是目录: {base_path}")
        # 获取目录内容并构建完整路径
    items = os.listdir(base_path)
    full_paths = [os.path.join(base_path, item) for item in items]
    return full_paths

def get_las_files(base_path):
    first_files = get_file_paths(base_path+'/first_file')
    second_files = get_file_paths(base_path+'/second_file')
    rgb_images_path = base_path + '/tif/rgbimg.tif'
    return first_files, second_files, rgb_images_path




def denoise_las(input_path, nb_neighbors, std_ratio):
    """
    用于对激光雷达点云数据（LAS格式）进行去噪处理，主要作用是移除点云中的离群点（噪声）
    :param input_path:
    :param nb_neighbors:
    :param std_ratio:
    :return:
    """
    las = laspy.read(input_path)
    points = np.vstack((las.x, las.y, las.z)).transpose()
    print(f'已加载 {input_path}')

    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(points)

    # 剔除孤立点
    cl, ind = pcd.remove_statistical_outlier(nb_neighbors=nb_neighbors,
                                             std_ratio=std_ratio)
    filtered = pcd.select_by_index(ind)

    print(f'已剔除孤立点 {input_path}')

    return np.asarray(filtered.points)

def generate_dem(point_arrays, res, dem_path):
    """
    从一组点云数据生成数字高程模型（DEM），并将其保存为GeoTIFF文件
    :param point_arrays:
    :param res:
    :param dem_path:
    :return:
    """
    # 1. 合并所有点云数据
    pts = np.vstack(point_arrays)
    # 2. 分离坐标分量
    x = pts[:, 0]
    y = pts[:, 1]
    z = pts[:, 2]
    
    # 3.计算点云边界范围
    xmin, xmax = x.min(), x.max()
    ymin, ymax = y.min(), y.max()

    # 4. 计算网格尺寸
    nx = int(np.ceil((xmax - xmin) / res))
    ny = int(np.ceil((ymax - ymin) / res))

    # 5. 创建规则网格
    xi = np.linspace(xmin, xmax, nx)
    yi = np.linspace(ymin, ymax, ny)
    xi, yi = np.meshgrid(xi, yi)

    # 6. 高程插值（核心步骤）
    from scipy.interpolate import griddata
    zi = griddata((x, y), z, (xi, yi), method='nearest') # nearest精度低，会出现阶梯状假块面，可改用linear但会很慢

    # 写入 GeoTIFF
    transform = from_origin(xmin, ymax, res, res)
    with rasterio.open(
            dem_path, 'w',
            driver='GTiff',
            height=zi.shape[0],
            width=zi.shape[1],
            count=1,
            dtype=zi.dtype,
            crs='+proj=utm +zone=48 +datum=WGS84',  # 坐标参考系，可根据实际需要修改
            transform=transform,
    ) as dst:
        dst.write(zi, 1)

    return zi


def get_file_point(files,tif_file_path):
    logging.info(f"get_file_point 开始")
    points_files = []
    for fp in files:
        pts = denoise_las(fp, denose_nb_neighbors, denose_std_ratio)
        points_files.append(pts)
    dem = generate_dem(points_files, resolution, tif_file_path)
    logging.info(f"get_file_point 结束")
    return dem




import numpy as np

def align_and_diff(dem1, dem2,
                   row_mode='pad_last',    # 'pad_first' | 'pad_last' | 'trim_first' | 'trim_last'
                   col_mode='pad_last',    # 同上
                   fill_value=np.nan,
                   threshold=-0.2):

    h1, w1 = dem1.shape
    h2, w2 = dem2.shape

    H = max(h1, h2)
    W = max(w1, w2)

    def adjust(arr, H, W, row_mode, col_mode):
        a = arr.copy()
        h, w = a.shape

        dr = H - h
        if dr > 0:
            # 需要填充
            if row_mode == 'pad_first':
                pad = ((dr, 0), (0, 0))
            else:  # pad_last
                pad = ((0, dr), (0, 0))
            a = np.pad(a, pad_width=pad, mode='constant', constant_values=fill_value)

        elif dr < 0:
            # 需要裁剪
            dr = -dr
            if row_mode == 'trim_first':
                a = a[dr:, :]
            else:  # trim_last
                a = a[:-dr, :]

        dc = W - w
        if dc > 0:
            # 填充
            if col_mode == 'pad_first':
                pad = ((0, 0), (dc, 0))
            else:  # pad_last
                pad = ((0, 0), (0, dc))
            a = np.pad(a, pad_width=pad, mode='constant', constant_values=fill_value)

        elif dc < 0:
            # 裁剪
            dc = -dc
            if col_mode == 'trim_first':
                a = a[:, dc:]
            else:  # trim_last
                a = a[:, :-dc]

        return a

    dem1_adj = adjust(dem1, H, W, row_mode, col_mode)
    dem2_adj = adjust(dem2, H, W, row_mode, col_mode)

    diff = dem2_adj - dem1_adj
    mask = diff < threshold # 只要下陷超过threshold

    return diff, mask



def scale_geotiff_to_mask_height(tif_path, target_height):
    with rasterio.open(tif_path) as src:
        arr = src.read()         # (bands, h0, w0)
        transform = src.transform

    img = reshape_as_image(arr)        # (h0, w0, bands)
    img_pil = Image.fromarray(img).convert("RGB")
    h0, w0 = img_pil.size[1], img_pil.size[0]

    scale = target_height / h0
    target_width = int(w0 * scale)

    img_resized_pil = img_pil.resize((target_width, target_height), Image.BILINEAR)
    resized_img = np.array(img_resized_pil)   # (H, W, 3)

    rows = np.arange(target_height)
    cols = np.arange(target_width)
    orig_rows = rows / scale
    orig_cols = cols / scale
    gr, gc = np.meshgrid(orig_rows, orig_cols, indexing='ij')  # 两个 (H, W) 数组

    rows_flat = gr.flatten()
    cols_flat = gc.flatten()
    xs_flat, ys_flat = xy(transform, rows_flat, cols_flat)
    xs = np.array(xs_flat).reshape(gr.shape)
    ys = np.array(ys_flat).reshape(gr.shape)

    # 构建 'lon,lat' 字符串形式的 coord_map
    H, W = gr.shape
    coord_map = np.empty((H, W), dtype=object)
    for i in range(H):
        for j in range(W):
            # coord_map[i, j] = f"{xs[i, j]:.6f},{ys[i, j]:.6f}"
            coord_map[i, j] = f"{xs[i, j]:.6f},{ys[i, j]:.6f}"

    # 提取黑色掩码
    white_mask = np.all(resized_img == 0, axis=-1)

    return resized_img, coord_map, white_mask




def overlay_mask_on_image(resized_img, mask, white_mask, overlay_color=(255, 0, 0), alpha=128):
    mask_flipped = mask[::-1, :]
    overlay_region = mask_flipped & (~white_mask)

    base_pil = Image.fromarray(resized_img).convert("RGBA")
    overlay = Image.new("RGBA", base_pil.size, (0, 0, 0, 0))
    draw = ImageDraw.Draw(overlay)

    ys, xs = np.where(overlay_region)
    for x, y in zip(xs, ys):
        draw.point((x, y), fill=overlay_color + (alpha,))

    composite_pil = Image.alpha_composite(base_pil, overlay)
    return composite_pil



def align_only(dem1, dem2,
               row_mode='pad_last',  # 'pad_first' | 'pad_last' | 'trim_first' | 'trim_last'
               col_mode='pad_last',  # 同上
               fill_value=np.nan,
               threshold=0.2):
    h1, w1 = dem1.shape
    h2, w2 = dem2.shape

    H = max(h1, h2)
    W = max(w1, w2)

    def adjust(arr, H, W, row_mode, col_mode):
        a = arr.copy()
        h, w = a.shape

        dr = H - h
        if dr > 0:
            # 需要填充
            if row_mode == 'pad_first':
                pad = ((dr, 0), (0, 0))
            else:  # pad_last
                pad = ((0, dr), (0, 0))
            a = np.pad(a, pad_width=pad, mode='constant', constant_values=fill_value)

        elif dr < 0:
            # 需要裁剪
            dr = -dr
            if row_mode == 'trim_first':
                a = a[dr:, :]
            else:  # trim_last
                a = a[:-dr, :]

        dc = W - w
        if dc > 0:
            # 填充
            if col_mode == 'pad_first':
                pad = ((0, 0), (dc, 0))
            else:  # pad_last
                pad = ((0, 0), (0, dc))
            a = np.pad(a, pad_width=pad, mode='constant', constant_values=fill_value)

        elif dc < 0:
            # 裁剪
            dc = -dc
            if col_mode == 'trim_first':
                a = a[:, dc:]
            else:  # trim_last
                a = a[:, :-dc]

        return a

    dem1_adj = adjust(dem1, H, W, row_mode, col_mode)
    dem2_adj = adjust(dem2, H, W, row_mode, col_mode)

    return dem1_adj, dem2_adj


def mark_color_blocks_with_number(pil_image,
                                  coord_map,
                                  target_rgb=(0, 0, 255),
                                  tol=10,
                                  min_area=100,
                                  max_area=10000,
                                  outline_color=(255, 0, 0),
                                  outline_width=2,
                                  text_color=(255, 0, 0),
                                  font_size=30,
                                  center_dist_thresh=10,
                                  font_path="DejaVuSans.ttf"):
    # 原图像尺寸
    W, H = pil_image.size

    img = np.array(pil_image.convert("RGB"))
    lower = np.clip(np.array(target_rgb) - tol, 0, 255).astype(np.uint8)
    upper = np.clip(np.array(target_rgb) + tol, 0, 255).astype(np.uint8)
    mask = cv2.inRange(img, lower, upper)

    num_labels, labels, stats, _ = cv2.connectedComponentsWithStats(mask, connectivity=8)
    draw = ImageDraw.Draw(pil_image)

    # 预加载固定字号字体
    try:
        font = ImageFont.truetype(font_path, font_size)
    except IOError:
        font = ImageFont.load_default()

    centers = []  # 存放最终 (idx, lon, lat)
    center_pixels = []  # 存放已绘制框的中心像素 (cx, cy)
    label_idx = 1

    for i in range(1, num_labels):
        x, y, w, h, area = stats[i]
        if area < min_area or area > max_area:
            continue

        cx = x + w // 2
        cy = y + h // 2

        skip = False
        for (px, py) in center_pixels:
            if np.hypot(cx - px, cy - py) < center_dist_thresh:
                skip = True
                break
        if skip:
            continue

        # 计算比例坐标
        rx1 = float(x / W)
        ry1 = float(y / H)
        rx2 = float((x + w) / W)
        ry2 = float((y + h) / H)

        # 绘制矩形框
        draw.rectangle([(x, y), (x + w, y + h)],
                       outline=outline_color,
                       width=outline_width)

        text = str(label_idx)
        bbox = font.getbbox(text)
        text_w = bbox[2] - bbox[0]
        text_h = bbox[3] - bbox[1]
        text_pos = (cx - text_w // 2, cy - text_h // 2)
        draw.text(text_pos, text, fill=text_color, font=font)

        lon_str, lat_str = coord_map[cy, cx].split(',')

        # 提取地理坐标
        lon, lat = float(lon_str), float(lat_str)

        # 记录信息到 dict
        centers.append({
            'index': label_idx,
            'center': (lon, lat),
            'topleft': (rx1, ry1),
            'bottomright': (rx2, ry2)
        })

        center_pixels.append((cx, cy))

        label_idx += 1

    return pil_image, centers

# 查看结果
# marked_img2.show()
# print("检测到的区域中心点：")
# for idx, lon, lat in center_list:
#     print(f"{idx}: 经度={lon}, 纬度={lat}")
#
# marked_img2.save('./static/marked_image2.png')
#
# print('地图的左上角坐标和右下角坐标', json.dumps([coord_map[0,0], coord_map[-1,-1]]))
#
# print('每个下塌框的序号和中心坐标', json.dumps(center_list))